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Abstract: The experimental bound on lifetime of U3, the neutrino mass eigenstate with 
the smallest component, is much weaker than those of ui and U 2 by many orders of 
magnitude to which the astrophysical constraints apply. We argue that the future reactor 
neutrino oscillation experiments with medium-baseline (~ 50 km), such as JUNO or RENO- 
50, has the best chance of placing the most stringent constraint on 1/3 lifetime among all 
neutrino experiments which utilize the artificial source neutrinos. Assuming decay into 
invisible states, we show by a detailed analysis that the U 3 lifetime divided by its mass, 
Ta/rus, can be constrained to be r^/ms > 7.5 (5.5) x 10“^^ s/eV at 95% (99%) C.L. by 
100 kt-years exposure by JUNO. It may be further improved to the level comparable to 
the atmospheric neutrino bound by its longer run. We also discuss to what extent U 3 decay 
affects mass-ordering determination and precision measurements of the mixing parameters. 
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1 Introduction 

Investigations on the possibility of neutrino decay has a long history, see e.g., [1, 2 ]. Since 
neutrino radiative decay is so tightly constrained [3], decay into invisible final states are 
more commonly discussed, for example, in the context of majoron models [4, 5].^ The 
bound on neutrino lifetime r depend on ( 1 ) whether daughter neutrinos are active or 
sterile [7], and in the former case (2) which neutrino mass ordering, normal or inverted, 
is realized. It also depends on (3) whether the neutrinos are Dirac or Majorana particles. 
However, in any one of these cases whenever the bound exists, its order of magnitude is 
given by the condition r/m ~ L/E for neutrinos with mass m and energy E that traverses 
distance L [ 8 ]. We refer the condition as the kinematic estimate. It is nothing but stating 
that decay effect is sizeable when traveling time t is comparable to lifetime riab of neutrinos 
in the laboratory (i.e., observer’s) frame. 

^ To be testable by various means we mention here, one has to arrange the majoron models such that 
neutrino lifetime is not too long [6]. 
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Astrophysical neutrinos, because of their long path lengths, are the promising sources 
for yielding the stringent bounds on neutrino lifetime. The solar neutrinos have been used 
to place bound on the lifetime of 1 ^ 2 , which is the dominant component of neutrinos 
under the assumption that i >2 decays into ui {vi) or into sterile states [8-12], In most 
cases the daughter neutrinos were assumed to be unobserved even if the decay into active 
neutrino is considered. The obtained bound is of the order of the kinematic estimate, 
7-2 7^2 10 s/eV. We note that the authors of ref. [11] argued that it is possible to 
constrain also the lifetime of vi by considering low energy pp and ^Be solar neutrinos, 
which have large ui components. Because of lower neutrino energies, the obtained bound 
for vi is better by about a factor of five than that for V 2 [Hj- 

Supernova neutrinos are potentially the most powerful source for bound on neutrino 
decay [13] which could lead to the bound r/m ~ 10^ s/eV according to the kinematic 
estimate. However, the available data is currently limited to the one which came from 
SN1987A [14, 15]. Moreover, the bound applies only to vi and/or V 2 which have large 
De components. Astrophysical neutrinos which have been observed by IceCube [16], in 
principle, are equally (or more) powerful as supernova neutrinos. But, to place the bound 
on r/m we need either identification of the sources, or complete determination of the 
neutrino flavor ratios [17]. Its typical order is estimated as r/m ~ 10^ s/eV assuming 1 
TeV neutrinos from AGN at 100 Mpc [ 8 ]. See [18] for recent discussions on the hypothesis 
of neutrino decays over cosmological distances. 

Leptonic decay of mesons leads to a bound on the majoron coupling constant with 
neutrinos gap (a,/3 = e,p,T) [19] of the order of ~ 10“^ — 10“^. For a comprehensive 
treatment of the majoron coupling bound with the pseudo-scalar as well as the scalar 
couplings, see e.g., [20] and an update [21]. When translated into 1^2 lifetime assuming 
decay into active ui, it gives a bound on r 2 /r?T ,2 in a range comparable with to somewhat 
weaker than the solar neutrino bound [ 8 ]. Using lepton decay channels the authors of [20] 
also obtained the bound of couplings which include Yla=e n r < 0 . 1 , which would 

lead to a less stringent bound on Ts/ms for active vs decay. 

Therefore, so far no stringent bound on 12 ^ lifetime appears to exist either from as¬ 
trophysical or from laboratory neutrinos. Probably, the best way to constrain 12 ^ lifetime 
would be to use neutrino oscillation phenomenon. It is because in this case one can se¬ 
lect the channels or region of kinematical phases that are sensitive to 123 decay effect. 
Then, it is natural to use the oscillation driven by Am 32 — hereafter referred to 

as atmospheric-scale neutrino oscillation for simplicity. We will see that despite quantum 
mechanical nature of the phenomenon the kinematic estimate applies. The bound on 2/3 
lifetime was obtained [23] by using the data collected by the Super-Kamiokande (SK) at¬ 
mospheric neutrino observation [24] as well as by the long-baseline accelerator experiments. 
See [25] for a recent update of the accelerator bound. 

In this paper, under the assumption of decay into invisible daughters, we analyze the 
bound on 1/3 lifetime which will be placed by the future medium-baseline reactor neutrino 
experiments such as JUNO [26] or RENO-50 [27] via observing the distortion of neutrino 
energy spectrum. In fact, we argue that they could provide, in principle, the most stringent 
bound on 12 ^ lifetime among all the artificial source neutrino experiments. It is because 
they can observe the effect of atmospheric-scale oscillations at the baseline around the 
maximum of the solar-scale oscillation, Tom = dvrU/Amli, which is longer than Lqm 
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for the atmospheric-scale oscillation by a factor of ~ 30. Assuming five years operation of 
JUNO, we obtain the bound ra/ms > 7.5 (5.5) x 10“^^ s/eV at 95% (99%) CL. See section 2 
for details, and for the relationship between our argument and the bound obtainable by 
using atmospheric neutrino data. 

The principal objective of the medium-baseline reactor neutrino experiments is to 
determine the neutrino mass ordering. Then, the immediate question is whether it would 
be disturbed if possibility of neutrino decay is taken into account in fitting the data, or 
more drastically, when 1/3 would actually decay. It was also noticed that such experiment 
has a potential of determining the mixing parameters in a high precision [28, 29]. In fact, 
the recent works with much more elaborate treatment of experimental errors reported an 
extreme precision of sub-percent level for sin^ 0 i 2 and Am^i [30, 31]. Then, the natural 
question is whether or to what extent these sensitivities could be affected when possibility 
of 1/3 decay is turned on. In section 6 we address these questions. 


2 Uniqueness of the medium-baseline reactor neutrino experiments 


In this section we try to convince the readers that JUNO/RENO-50 is, in principle, the 
highest sensitivity experiment among all those which utilize artificial source (or beam) 
neutrinos in detecting the possible decay effect of 1 ^ 3 . In this paper, we consider the case 
of invisible decay of U 3 . That is, we assume that the decay products are either some sterile 
states, or can involve active and/or 1^2 state but with significant energy degradation 
such that daughter neutrinos cannot be observed. The latter possibility necessitates the 
neutrino mass ordering to be the normal type. 

When i-th mass eigenstate neutrino decays with lifetime r* at rest, the energy Ei 
(more precisely, the energy difference normalized appropriately) of propagating i-th mass 
eigenstate neutrino can be written as 



where 


( 2 . 1 ) 


( 2 . 2 ) 


is a Lorentz dilated lifetime. 

As will be shown in appendix A the decay of U 3 produces the following two character¬ 
istic modifications in the oscillation probabilities in vacuum:^ 

• Reduction of the atmospheric-scale oscillation amplitude with the form cos I —® ^ • 

• Decrease of normalization of the probabilities by the amount proportional to (l — . 


Here, L is the distance traveled by neutrinos. To make the point clearer we have used an 
approximation ~ Am |2 = where Am^j = m| — mf. The first feature stems 

from the fact that the effect comes from the interference between the first—second mass 
eigenstates and the third. The second feature represents the effect of decaying z-th mass 
eigenstate projected onto the initial and final neutrino flavor states. It is independent of 

^ The matter effect does not appear to affect the conclusion of this section in a significant way. 
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neutrino oscillation and exists even in the limit of vanishing Am^j, as it must. This effect 
is prominent in (and i>^) disappearance channel, but is negligible for Ue (and Ve) survival 
probability. Whereas the first effect is equally important in the all oscillation channels, as 
will be shown in appendix A. 

The decay effect is negligible for baseline L with T^L <C 1, and it becomes significant 
only when TsL > 1. In most of the oscillation experiments which use man-made neutrino 
sources, the baseline is set to the first oscillation maximum of atmospheric-scale oscillation, 
Aml^L/i 2 E) ~ TT. In such setting, it is likely that we obtain the lower bound on the 
width Ta as 


1 ^ IttE 

Ts-^- Am2 


atm 


(2.3) 


From (2.2), we know that Tzjvaz = (^sE) Using this relation in (2.3) one can estimate 
the achievable upper bound on ra in such setting as 


T3 ^ 2tt /_ L\ 
ms ~ V Ej 


(2.4) 


where the last equality in parenthesis assumes that the distance L is at the first oscillation 
maximum. It implies that the kinematic estimate applies. Using (2.4), one can estimate 
the upper bound on rs/ma at 


^ ~ 1.7 X 10-12 f -^ 

ma V 2.4 X 10“3 eV^ 



(2.5) 


In the JUNO or RENO-50 setting, hereafter referred simply as JUNO setting for sim¬ 
plicity, the detector is placed approximately at the maximum of the solar-scale oscillation, 
Aml^L/{2E) ~ TT. Then, we can achieve about 30 times longer lifetime bound 


— ~ 5.5 X 10-11 
m 3 


Am^ 


21 


7.5 X 10-5 eV2 


s 

eV' 


( 2 . 6 ) 


It should be stressed that the JUNO setting is unique (among artificial neutrino source 
experiments) in making observation of the atmospheric-scale neutrino oscillation at the 
distance of the solar-scale oscillation maximum possible. This completes our argument that 
JUNO is the highest sensitivity experiment among (practically) all the ongoing or proposed 
artificial-source neutrino experiments as far as the 1^3 lifetime bound is concerned. 

Of course, the above argument does not necessarily imply that the JUNO bound on 1/3 
lifetime must be the severest one achievable by all the neutrino experiments. The likely (and 
probably unique) exception is the one placed by observation of the atmospheric neutrinos. 
The naive kinematic estimate for the 1^3 lifetime would lead to the sensitivity to 


T 3 L 
_ _ 

m 3 E 


10^ km 
1 GeV 


~ 3.3 X 


s 


(2.7) 


which suggests that the bound by the atmospheric neutrinos is comparable as the one 
by JUNO. However, in the case of atmospheric neutrino experiments they can observe 
neutrinos in much wider energy and baseline ranges than those used in the above estimate. 
Therefore, we expect that the bound on ra/ms from the atmospheric neutrinos is tighter 
than the naive estimate given in (2.7). In fact, the authors of ref. [23], by using the SK 
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atmospheric neutrino data and the others, obtained the bound T-ijraj, > 9.3 x 10 s/eV 
at 99% CL, which is stronger than the one shown in Eq. (2.7) by about a factor of 3.^ 


3 Effect of z /3 decay on the oscillation probabilities and the observable 


The i^e survival probability relevant for reactor neutrinos for the baseline L in vacuum is 
given, under the approximation ~ Am |2 = 


^atm5 


as 


P{ve —^ = 1 — Ci3 sin^ 2012 sin^ 


AE 


_ ^4 


(l_e-r3L) - Isin^ 20i3 


^ / Am^i^L \ rgL 


(3.1) 


where Cij = cos Oij and Sij = sin Oij. See appendix A for derivation. For simplicity and as 
a good approximation, we ignore the matter effect in this work. 

As we stated in the previous section, there are two types of terms which are affected 
by the neutrino decay. However, they come with vastly different magnitudes in the T'e 
channel; since the current neutrino data implies ~ 0 . 02 , the coefficients of (l — 
and cos terms are, respectively, S 13 ~ 5 x 10“^ and ^ sin^ 20 i 3 ~ 4 x 10“^. 

Therefore, the former oscillation-independent Ve attenuation effect should be negligible. 

Then, the question is how such decay-affected oscillation probability manifests itself 
into the observable quantities. To give a feeling to the readers we show in the upper panel 
of Fig. 1 the energy spectrum of events as a function of positron deposited energy, or visible 
energy Fvis calculated by convoluting the neutrino flux and the cross section of the inverse 
/?-decay (IBD) + n) reaction. The three curves with different combinations of 

colors and line types correspond to the following three cases: no decay (black solid line), 
Tzjni-i = 5 X 10“^^ s/eV (red solid line), and T-ijm-i = 10“^^ s/eV (green dashed line). See 
appendix B for the procedure to obtain these curves. As we can see from the upper panel 
in Fig. 1, the decay of 1^3 state tends to average out the fast oscillation (wiggles) driven by 

Since decreasing 0 i 3 and 2/3 decay both reduce the amplitude of atmospheric-scale 
oscillation, they could potentially be confused with each other. Fortunately, the confusion 
is precluded by the precision measurement of 0 i 3 done by the short baseline reactor neutrino 
experiments. (See the related comment at the end of section 5.) Currently, sin^ 20 i 3 is 
measured with the uncertainty of ~ 6% [34], while a possibility of reaching the ultimate 
error of ~ 3% by the end of 2017 is mentioned in [35]. 

Degradation of the oscillation amplitude can also occur by a finite energy resolution 
of the detector, which would cause another confusion with the decay effect. To illustrate 
this point, we show in the lower panel of Fig. 1, the difference between the energy spectra 
in the absence and in the presence of vj, decay for two cases of energy resolution, 3% (red 

® For comparison we note that the kinematic estimate for JUNO, (2.6), is different from our results 
based on analysis by ~ 40% only. It should also be noticed that if a re-analysis in [23] would be done 
with the currently accumulated SK data, the bound should become even tighter. It is also more stringent 
by a factor of 1.7 than our JUNO five-years bound to be obtained in section 5. 

^ The visible energy ifvia is approximately related to neutrino energy E as ifvis — E — (m„ — rrip) -f me, 
where m„, nip, me, are, respectively the mass of neutron, proton and electron. 

® In this connection we note that quantum decoherence [32] has the similar effect on energy spectrum 
of Ue, and the relevant phenomenology is discussed with the JUNO setting in ref. [33]. 
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Figure 1. Upper panel shows the expected event distribution as a function of the visible energy 
of positron at the 20 kt detector placed at L = 52.5 km from the reactor complex of 35.8 GW 
thermal power. The case without decay effect is indicated by the black solid curve whereas the case 
with decay effect is shown by the red solid and green dashed curves, corresponding, respectively, 
to T3/TO3 = 5 X 10“^^ and 10“^^ s/eV. We also show the individual contributions coming from 
the reactor complex at the medium-baseline L = 52.5 km (blue dotted curve), from the far reactor 
complexes located at Daya Bay with L=215 km and Huizhou with L = 265 km (violet dashed 
curve) and geoneutrinos (orange dash-dotted curve). In the lower panel, in order to see clearly the 
importance of the energy resolution, we show, for T^jmz = 5 x 10“^^ s/eV, the difference of the 
cases without and with decay effect, or dN/dE (without decay) - dN/dE (with decay), as a function 
of E^is for 3% and 6 % energy resolution, respectively, by the red solid and blue dashed curves. 

solid curve) and 6 % (blue dashed curve) at 1 MeV. Roughly speaking, doubling the energy 
resolution causes an amplitude attenuation of the signal by a factor of three. We observe 
that the detectability of the decay effect is strongly dependent on the energy resolution. 
Fortunately, the 3% energy resolution is to be reached by JUNO and RENO-50 [26, 27]. 

To our understanding, the remaining questions which need to be addressed are as 
follows: 

• As shown in Fig. 1, the decay effect is to reduce the atmospheric-scale oscillation 
amplitude. This is the oscillation to be utilized to determine the mass ordering in 
JUNO. Then, the obvious question is to what extent the mass ordering determination 
could be disturbed if the effect of t's decay is taken into account. 

• Another relevant question would be to what extent the sensitivities to mixing param¬ 
eter measurement in JUNO could be disturbed by 1^3 decay. 

We will discuss these issues in section 6 from the following two different viewpoints (as¬ 
sumptions): (1) There is no decay in the input data set, but we consider the decay effect in 
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the output (fit), and (2) 1/3 actually decays with the lifetime which is marginally consistent 
with the current and the expected JUNO bound. 


4 Analysis method 

For definiteness, throughout the paper (including Fig. 1), we assume that the true (input) 
values of the oscillation parameters are given as follows, 


Amli = 7.50 X 10“®eV^ sin^ Ou = 0.304, 
= 2.46 X 10"^eV^ sin^ (9i3 = 0.0218, 


(4.1) 


which are taken from the best fitted values of the one of the recent global analysis [36], 
and assume the normal mass ordering (Am|^ > 0), nnless otherwise stated. 

In our statistical analysis, we formally divide the fnnction in three terms as 


2 _ 2 I 2 

X — Xstat 4” Xr 


+ Xs 


(4.2) 


^param ' Asys’ 

as done in [30, 31]. The first term, compnted by taking the limit of infinite number 

of bins which is justified because of the large number of events expected at 20 kt detector, 
~ 1.4 X 10^ for 5 years of operation, as follows [30, 31], 

2 


/dN' 


obs 


Xstat 


rf 


dEv 


dE^ 


Y. 


dNf^\ 


z=reac, U, Th 




/ C^ATobs 


V 




(4.3) 


/ 


where dN°^^/dEvis is the event distributions of the observed (simnlated) signal, and is 
the flux normalization parameters for reactor neutrinos as well as for geonentrinos to be 
varied freely snbject to the pnll term in (see below) and we integrate np to E^^^ = 8 
MeV. In onr analysis, following [31], we include the contributions not only from Yangjiang 
and Taishan reactor complexes at L = 52.5 km (approximated by a single reactor with the 
thermal power of 35.8 GW) but also the contributions coming from the reactor complexes 
located at Daya Bay {L = 215 km) and Huizhon (L = 265 km), as well as geoneutrinos. 
See appendix B for details. 

The second term in (4.2) takes into acconnt the current uncertainties of the standard 
oscillation parameters and is given by 


Xparam 


4 

E 

2 = 1 


fit 


Xj — X. 


a{xi) 


(4.4) 


where Xi and xf* (i = 1-4) denote, respectively, the assumed true (input) and fitted values 
with xi = sin^di 2 , X 2 = X 3 = sin^dis, X 4 = Am|;^. For the values of a{xi), we 

take the current 1 sigma nncertainties determined by the global fit [36], (T(sin^ Ou) = 4.1%, 
a(Am 2 i) = 2.4%, cj(sin^di 3 ) = 4.6% and cr(Am|;^) = 1.9%. 

The last term in (4.2), Xgyg, takes into account the contributions of two kind of exper¬ 
imental systematic uncertainties we consider 


X^ = 

Asys — 


fit 

reac 






cht \ ' 

— ) 


-k 




fit 

Th 




-k 


fit 


'n_ 

(Jn 


(4.5) 
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where and describe, respectively, the normalization uncertainties for reactor 

neutrinos, uranium and thorium induced geoneutrinos. For them we take = 3% [37, 
38] for reactor neutrinos and = 20% for geoneutrinos following ref. [31]. In 

addition to the normalization uncertainties, we also consider the uncertainty of the energy 
resolution by including the last term in (4.5) with arj (see below). 

For simplicity, we assume 100% detection efficiency, Cdet = 1- The effect of uncertainty 
in the detection efficiency (IBD selection, fiducial volume, etc) is thus neglected, since it 
is difficult to reliably estimate. However, those errors, being energy independent, should 
hardly affect the sensitivity to shape-modulating effects caused by the neutrino decay as 
well as by the differing mass orderings. 

With regard to the energy resolution, we only consider the stochastic term, i.e., cte/E = 
0.03(1 -|- T])/E/M.eY contribution where r] is the parameter accounting for the energy 
resolution variation.® The non-stochastic term(s) are expected to be the most relevant, even 
dominant, at energies > 3 MeV. However, this is a complex experimental matter impossible 
to anticipate at this stage, so its impact is simply neglected here. Some deterioration of the 
sensitivity should be expected, if those terms were considered. Because of the illustrated 
dependence of the energy resolution on the decay sensitivity, its uncertainty is taken into 
account by using a pull term with = 10%.^ We feel it appropriate to consider no other 
systematic error on the energy scale at this stage. 

The is computed in the following way: In order to derive the expected bound on the 
1/3 decay lifetime, for our input data, we consider the JUNO setup to compute /dE^is 

assuming no decay {rs/ms = 00 ) using the oscillation parameters given in (4.1). Then we 
try to fit such dN°^^/dEvis by minimizing the function varying freely 9 i2, Arn^i, ^ 13 , 
^reac) Cuj ^Th; d, and Tslm^ subject to the pull terms in (4.4) and (4.5). In this 
way we can estimate the sensitivity to Tsjm^ and derive the bound on it, and at the same 
time, can determine the allowed regions of other parameters as well. 

Using the function, we will determine the allowed range of Tsfms by the condition 

= X^ - xLin = 2.71, 3.84 and 6.63 (1, 4 and 9), (4.6) 

at 90%, 95% and 99% (1, 2 and 3 a) CL, respectively, for one degree of freedom. For the 
case where we show the allowed regions in the plane spanned by any combination of two 
out of nine parameters, we use the condition Ax^ = 2.3,6.18 and 11.83, respectively, for 1, 
2 and 3 a CL for two degree of freedom. We note that Xmin = 0 by construction because 
we do not take into account the statistical fluctuation in simulating the artificial data. 

5 The bound on T^/m^ 

In this section, we derive the bound on T-^jraj, assuming that 1^3 decays into invisible states, 
whereas ei and 1^2 are stable. For definiteness, we assume the normal mass ordering. But, 
given the expression of P{ue in (3.1), we expect that our results are valid also for 

the case of inverted mass ordering. 

® See appendix B for the Gaussian energy resolution implementation in our analysis. 

^ Notice that the energy resolution uncertainty, usually estimated by using calibration source data, 
has never been a sensitive quantity for physics parameter determination so far, but it must be carefully 
considered by experimental collaborations in completion of the detector design. 
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As our standard setup, we assume 5 years running of the JUNO detector with the 
fiducial volume of 20 kt which is placed at L = 52.5 km away from an effective single 
reactor with the thermal power of 35.8 GW. We assume the detector’s running with 100% 
efficiency, Cdet = 1- But, if Cdet < 1 the running time must be scaled to [5/edet] years to 
obtain the same results. Since JUNO may take data for a longer period, we also consider 
the case of exposure for 15 years. See appendix B for details of our calculation. 



Figure 2. — Xmin B shown, by the red (blue) curves for 5 (15) years of data taking, 

as a function of the fitted value of rs/ma calculated for the JUNO detector placed at L = 52.5 
km from a reactor with 35.8 GW thermal power, assuming 5 years of exposure and 100% detection 
efficiency. We have taken that the true (input) value of T 3 /TO 3 is infinite (stable 1 / 3 ). The solid curves 
correspond to the results obtained by using our full defined in (4.2) whereas the dashed ones 
correspond to the case without assuming systematic errors. The contributions from the reactors 
at Daya Bay and Huizhou as well as those from geoneutrinos are taken into account. The bound 
comes from the SK atmospheric neutrinos plus long-baseline oscillation experiment obtained in [23] 
is also indicated by the vertical black dashed line. 

In Fig. 2 we show — Xmin ^ function of the fitted value of Ts/m^ where 

the input (true) value of T^frus is assumed to be infinity, i.e., is stable, with oscillation 
parameters given in (4.1). All the other eight parameters are marginalized in the fit. The 
case of exposure for 5 (15) years is shown by the blue (red) solid curve. To exhibit the 
effect of the systematic error onto the bound, we also show by the dashed curves Ax^ for 
the case without all the systematic errors. 

From our result shown by the solid blue curve in Fig. 2, we can conclude that if the 
data is consistent with no-decay hypothesis, the range 

— < 7.5 (5.5) X 10“^^ 47 (5-1) 

m 3 eV 

can be excluded at 95 (99%) % CL by 5 years exposure by JUNO. 

As mentioned in section 2 , this bound is about a factor of 1.7 weaker than the bound 
rs/ma > 9.3 x 10“^^ s/eV at 99% CL based on the data of atmospheric neutrinos (as well as 
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of long-baseline oscillation experiments) [23]. The bound is indicated by the vertical black 
dashed line in Fig. 2. By considering an extended running of 15 years, the JUNO bound 
can be improved to 11 (8.5) xlO“^^ s/eV at 95 (99%) % CL, which is barely comparable 
to the atmospheric neutrino bound. A comparison between the bound of the order of 
^ 10“^° s/eV obtained here and in [23] to the one deduced by using lepton decay 
channel [20] is described in [23]. They obtained the bound on 1/3 — majoron coupling 
ISsSp ^ 10“^ (?7is/l eV)~^ at 90% CL assuming m 3 3 > m^, where m^ denotes the sterile 
neutrino mass. 

We also considered the hypothetical situation where @13 were unknown at the time of 
JUNO running (not shown). We did it by removing the pull term for @13 from our 
and found that the bound on T^lmz would become worse by about a factor of three. It is 
because the 1^3 decay and the effect coming from the uncertainty in 0 i 3 can be confused 
with each other, as we mentioned in Sec. 3. 


6 Impact of z /3 decay on the determination of the mass ordering and the 
oscillation parameters 


6.1 Impact of t's decay on the mass ordering determination 

Now, we address the question of effect of 1^3 decay onto the mass ordering determination 
in JUNO/RENO-50, as promised at the end of section 3. We use the same analysis tool as 
used to obtain the bound on lifetime in the previous section. The input mass ordering 
is always taken to be normal in this section. 



Figure 3. — Xmin shown as a function of the fitted value of Amgg. The solid blue 

and red curve correspond, respectively, to the case where the fit is performed assuming the normal 
(right) and inverted (wrong) mass ordering (MO) for the standard oscillation. The blue (red) dotted 
and dashed curves corresponds, respectively, to the case where the input value of Ta/ma is 00 and 
10“^° s/eV for the case of normal (inverted) mass ordering. 


To know the effects of z /3 decay on the resolution capability of the mass ordering 
in JUNO, we compute Ay^ = y^ — Xmm taking both the normal (right) and in- 
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verted (wrong) mass ordering.® In Fig. 3, Ax^ is plotted as a function of ob¬ 

tained by marginalizing all the other parameters. For the abscissa in Fig. 3 we use 
Anigg = |c? 2 Am|i + sf 2 Am| 2 |, which we believe to be the appropriate variable to dis¬ 
cuss resolution of the mass ordering in medium baseline reactor experiments [39]. It is 
proposed as the effective atmospheric Am^ determined by Ue disappearance experiment in 
vacuum [40] , which agrees in a good approximation with the one measured by the reactor 
013 experiment [41]. 

The blue and the red curves are for the normal (input) and the inverted (wrong) mass 
orderings, respectively. The three line types, which are common to the both mass orderings, 
correspond respectively to: 

(i) Solid curve: no u decay in both the input and the fit, the case of standard oscillation, 

(ii) Dotted curve: no v decay in the input but 1/3 decay is allowed with prior-unconstrained 
lifetime in the fit, and 

(iii) Dashed curve: 1^3 decay with lifetime Tsims = 10“^*^ s/eV is assumed in the input and 
allowed in the fit. 

The global features of the results can be summarized as: 

• From (i) to (ii) there is only minor change (up to ~ 1) in Ay^. It means that allowing 
the decay only in the ht little affects the output. 

• From (i) to (iii) there are appreciable changes in Ay^ in a manner which depend very 
much on the mass ordering. In the case of normal (right) mass ordering the change 
in Ax^ upon turning on 1/3 decay both in the input and output is modest, slightly 
opening up the Gaussian parabola. 

• If Ax^ difference between the right (input normal) and the wrong mass orderings 
without decay is denoted as ^xio-decay^ difference at the minima across the 
different mass orderings becomes ~ AXno_dg(,ay — 5 when the 1^3 decay is turned on 
with lifetime comparable to the current and the JUNO bound, r^jms = 10“^^ s/eV. 

Thus, we have observed that the U 3 decay has a big impact on mass ordering resolution in 
JUNO, significantly worsening the sensitivity. 

The readers must be noticed that we carefully avoided to make a quantitative state¬ 
ment on how large is Ax„Q_decay> restricted ourselves into the change due to the decay 
effect. It is likely that our procedure to simulate the distortion of the event number energy 
distribution is insufficient to reliably extract the absolute confidence level for the mass 
ordering determination. In particular, we do not take into account the uncertainties re¬ 
lated to the non-linearity of the energy measurement, whose control would be the key to 
the success in the mass ordering measurement in JUNO.® On the other hand, the omitted 
error in computing x^ may not affect our discussion of the effect of decay in a significant 
way, because it affects the spectrum as simple reduction of the oscillation amplitude.^® 

® At the end of this section 6.1, we will place a long clarifying remark on how to interpret onr Ax^- 
® This issue was raised in [39], and the crucial importance of controlling the energy-scale errors on mass 
ordering determination was illuminated with experimental perspectives by the authors of ref. [42]. 

In addition, there is a snbtle issues related to statistical treatment of mass ordering determination. In 
the nsual case one can attribute to \dKxf~ the meaning of la significance. The special feature that the 
fitting parameter is a discrete variable requires a different treatment to evaluate correctly the CL for the 
mass ordering determination [31, 43, 44]. Since we do not intend to elaborate this point, and because of 
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6.2 Impact of decay on the determination of the oscillation parameters 

We briefly discuss in this section possible effects of t's decay on determination of the stan¬ 
dard oscillation parameters in JUNO as well as on the values of output uncertainties. 
Detailed features of parameter correlations with and without r's decay, including the cor¬ 
relations with the systematic errors, will be discussed in appendix C. As in the previous 
section 6 . 1 , we consider the three cases, (i) standard oscillation, (ii) decay effect only in 
the fit and (hi) decay effect both in the input and in the fit. For the finite input value of 
lifetime we always use ra/ma = 10 “^*^ s/eV. 

To present our results in a clear way, we always use the following line symbols through¬ 
out figures 4 (this section), 5 (appendix C.l), and 6 (appendix C.2): For the case (ii) we 
show the 1, 2 and 3cr allowed regions, respectively, by the filled blue, yellow and green 
colors. For the case (i) the allowed regions are delimited by the black solid curves, and for 
(hi) by the black dotted curves. 

In Fig. 4 we show the allowed regions of the parameters in (a) sin^ 612 - (left 

panel), and (b) sin^ 6*13 - (right panel) space. We observe in the left panel that there 

is no visible effect of the decay on the determination of the parameters in the 1-2 sector 
of the MNS matrix, sin^ 0 i 2 and Arn^i- Whereas in the right panel, it is revealed that 
accuracy in measurement of (of sin^ @13 also but only slightly) is affected when 

decays. 



Figure 4. The regions allowed at 1, 2, and 3cr CL are drawn, respectively, by the filled colors of 
blue, yellow and green in the parameter space spanned by (a) sin^ 612 - Am^i (left panel) and (b) 
sin^ 6 »i 3 - Am^i ( right panel) for the case (ii) where no decay is assumed in input but allowed in 
the fit. For comparison, the cases of (i) the standard oscillation fit without decay and (iii) the case 
with the input value of T 3 /TO 3 = 10“^° s/eV are shown, respectively, by the black solid and black 
dotted curves. 5 years of exposure is assumed. 


To obtain a more comprehensive view of effect of decay and to represent the change 
in the sensitivities in a more quantitative way, we present in Table 1 the fractional uncer¬ 
tainties (in %) of all the parameters determinable in the three cases (i-iii) in the JUNO 
setting. In doing so we do not restrict to the oscillation parameters, but also include the 
ones which describe the systematic errors. From Table 1 and Fig. 4 we observe: 

the crude nature of our construction, we recommend the readers to use the information presented in 
Fig. 3 only to have a feeling on how 1/3 decay affects the determination of the mass ordering in JUNO. 
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Table 1. Fractional errors in percent (at la) of the oscillation as well as the parameters for 
systematic uncertainties, before (prior) and after the fit to 5 years of JUNO data with 100% detection 
efficiency, for the case where the true mass ordering is normal. They are determined by the condition 
= 1 for one degree of freedom. We consider three cases (i) standard oscillation fit without 
decay (ii) decay sensitivity fit (with true T^jm^ = oo) and (iii) decay measurement fit (with true 
Tsims = 10“^° s/eV). 


parameter 

prior error (%) 

(i) 

fitted error (%) 

(ii) 

(iii) 

sin^ 012 

4.1 

0.35 

0.35 

0.35 

Amf2 

4.1 

0.21 

0.21 

0.21 

sin^ 013 

4.6 

3.7 

3.8 

4.3 

Amfs 

1.9 

0.12 

0.12 

0.16 

1 T Creac 

3.0 

0.50 

0.50 

0.51 

l+^u 

20 

12 

12 

12 

1 + ^Th 

20 

13 

13 

13 

l+r] 

10 

5.5 

6.0 

7.1 


• Comparison between the columns (i) and (ii) indicates that if the decay effect is 
considered only in the fit, its impact is virtually absent except for @13 and the energy 
resolution uncertainty parameter rj. 

• Comparison between the columns (i) and (iii) tells us that if the decay effect is 

considered for both in the input and the fit, there is appreciable impact of the size 
of approximately 10%-30% but only for sin^6*i3, and l+rj. 

As a whole, we conclude that the impact of the decay on the determination of the mass 
and mixing parameters is rather small. 

7 Conclusions 

In this paper, we have investigated the question of how strong a constraint on decay lifetime 
of the massive neutrino state r'3 can be placed by the medium-baseline (L ~ 50 km) reactor 
neutrino experiments, JUNO or RENO-50, which we referred simply as JUNO in most part 
in the text. Assuming decay into invisible states and ti,T 2 T3, we found that the bound 
Tslfns > 7.5 (5.5) x 10“^^ s/eV at 95% (99%) CL can be obtained by JUNO with its five 
years exposure with 100% efficiency. 

In fact, there is a simple reason why the JUNO setting can offer the chance of deriving 
the most stringent bound on the neutrino lifetime among all the experiments which utilize 
neutrinos from the artificial sources. Given that such experiments are designed to have 
sensitivities at around the first oscillation maximum (in exploring the atmospheric-scale 
neutrino oscillations), the kinematic estimate implies r/m ~ LjE ~ 2t:/ However, 
if there is such an experiment that can explore atmospheric-scale oscillations at the distance 
of solar-scale oscillation maximum, the bound on r/m would become severer by a factor 
of Am^tm/Arnl]^ ~ 30. This is what JUNO does. 

The bound on we obtained for 5 years exposure of JUNO has the same order 

of magnitude as (but is somewhat weaker than) the one obtained with the atmospheric 
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neutrino data. We have examined the question of whether the 15 years running of JUNO 
can tighten up the bound to the level of the current atmospheric neutrino bound, and 
obtained an affirmative answer. It is nice to see that the comparably strong lifetime 
bound can be deduced by the two completely different experiments, one observing reactor 
neutrinos at the medium baseline, and the other measuring the atmospheric neutrinos. 

We have also discussed to what extent the t'a decay could affect the determination of the 
mass ordering as well as the precision parameter measurement by JUNO. We considered the 
two cases (in the notation in section 6), (ii) the decay effect is absent in the input data set 
but allowed in the output (fit), and (iii) 1/3 decay in input with the lifetime rs/ma = 10“^° 
s/eV, which is marginally consistent with the current and our JUNO bound. For the 
mass ordering determination, we found the significant impact of the decay (reduction of 
Ax^ by five units) but only for the case (iii). With regard to influence of decay to the 
measurement of the oscillation parameters, we found that the impact is rather small except 
for Amg]^ again only for the case (iii), in which the uncertainty of Arugg would become ~ 
30% larger. 


A Neutrino oscillation probabilities in the presence of neutrino decay in 
vacuum 


The neutrino evolution in vacuum with neutrino decay is governed by the usual Schrodinger 
equation by replacing the neutrino energy Ei by the one in (2.1). The S matrix whose 
elements describe neutrino flavor transition in vacuum as Vaix) = <S'ay3(x)r'^(0) is given by 

Sa^ix) = Y, (A.l) 

i=l 


where x is the distance traveled by neutrino and the notations for Uai, Trii, E, and Tj are 
the same as in section 2. Then, the general expression of the three-flavor expression of 
oscillation probability with neutrino decay in vacuum is given by [45] 


P(z //3 Ua) = |5'«/3| = P(i //3 z/q, : no decay) - Y\Uaif\Upi\'^ (l “ e 


r^+r ■ 


2 Y COS Aji ( 1 - e 2 ' ^ j +‘^Y sin Ajt ( 1 - e 2 "^), (A.2) 

j>i ^ ^ j>i 


where ^ OLiUj^iU(y^jJJ ^ j ^ ^ji = Am‘jj^x/{2E). The first term in (A.2), the oscillation 

probability in the absence of neutrino decay, is given by the familiar expression, 

P{vp z/q, : no decay) = 5ap - 4 ^ Re[Ca,/ 3 y] sin^ (^ ) “ 2 ^ ImfCa/jjj] sinAjj. (A.3) 

j>i ^ 2 


If the approximations Tix ~ T 2 X <C 1, and Tax ~ 1 holds we can simply (A.2). The 
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disappearance oscillation probability takes the particularly simple form as 


= 1 - 4|[/„i|^|C/a2p sin^ 


i21 


4 |C/„i|Va3psin2 


i=l,2 


A 


3 i 


- lUasl^l - - 2 Y, \Uai\^\Ua3\^cosA3i (l - , 


i=l,2 


« 1 - 4\Uai\‘^\Ua2\^ Sin^ (^) ~ (1 - \Ua3\‘^) sin^ 

- |C/„3|^ (1 - _ 2|C/„3|2 (1 _ |?7^3|2) cos Aatm (l " e"^") , (A.4) 

where in the last line we have used the approximation A31 ~ A32 = Aatm-^^ In the 
appearance channels {a ^ f3), the oscillation probability also simplifies: 

P{u^ Va) = Pii^p : no decay) - 117«3p| 17/33(l - 

- 2 ^ Re[Cai3ij] cos A3i (^1 - ^ X] sin As* (^1 - 

i=l,2 i=l,2 


^ P{U 0 : no decay) - |C/a 3 |^|^ 7 / 33 p (l - e 

- 2|C/„3p|t//33pcosAatm(l-e"^) (A.5) 

Under the approximation A 31 A 32 = Aatm and because of anti-symmetry of Im[C/ojC/^j?7*j?7^j] 

(= lTa.[Cai 3 ij]) under i -H- j, CP violating term in the decay-width dependent part of the 
oscillation probabilities vanishes. 

Some remarks are in order: 


• In disappearance channel the coefficients of (l — and cos Aatm terms are 

given by \U^ 2 ,\^ = ~ 0.24 and 2\U^z\'^ (l — |U/j 3 p) — 0.50, respectively. There¬ 

fore, the oscillation independent depletion effect may be as important as the effect 
of diminishing amplitude in the atmospheric-scale oscillations. 

• In Vf, appearance channel the corresponding coefficients of the two r 3 -dependent 
terms are given by \Ue 3 \‘^\U^z\‘^ = S 23 '^i 3 Si 3 — 1 - 1 x 10 “^ and twice of that ~ 2 . 2 x 10 “^, 
respectively. They are comparable with each other, and are similar in magnitude as 
the main oscillation term in vacuum. 


B Event energy spectrnm 


The distribution of the number of events coming from the inverse /3-decay (IBD) reaction, 
Ue + P ^ + n, as a function of the visible energy is given by, 


diV(^vis) 


— n^ptexp 


Ee / dE 

'J Emi-, 


E 


dcl)i{E) da{E^,Ee) 

^deti-^e 


2=reac,geo—z/ 
^Piiye ^ ^ei Li") 1 


dE 


dEf, 


(B.l) 


If we were to go through a theoretical discussion of the effect of neutrino decay on determination of 
the neutrino mass ordering, we should have made our treatment more elaborate one. It would necessitate 
use of the effective atmospheric ArUee as mentioned earlier. However, determination of mass ordering is 
little affected by the decay effect as far as rs/ms is much smaller than the current bound. Though there 
is some influence if the true value of ra/ma is comparable with the current bound, we prefer to remain in 
the simpler treatment above in this paper. If necessary, it is straightforward to formulate this problem by 
taking the framework adopted in [39]. 
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where Up is the number of target (free protons), texp is the exposure, edet is the detec¬ 
tion efficiency, d(l)i{E)/dE is the differential ffuxes of reactor neutrinos and geoneutrions, 
da{E^, Eg)/dEg is the IBD cross section, Pi{i'e Ue',Li,E) is the i^e survival probabilities 
and R{Ee, E^is) is the Gaussian resolution function (see below). 

For simplicity, we set Cdet = Ij and as a reasonable approximation for our purpose, we 
ignore the neutron recoil in the IBD reaction and simply assume that neutrino energy, E, 
and the positron energy, E^, is related as E^ = E — {rrin — rup) ~ FI — 1.3 MeV. Due to 
the finite energy resolution, the event distribution can not be obtained as a function of Eg 
(true positron energy) but as a function of the reconstructed or so called visible energy, 
Evis, which is approximately related to neutrino energy as Svis — E — {mn — rUp) -|- rUg, 
after taking into account the energy resolution (see the text below). Regarding the cross 
section, da{Ei,,Ee)/dEe, we use the one found in [46]. 

For the JUNO detector, we assume the same proton fraction ~ 11% as the Daya Bay 
detectors [47] which implies ~ 1.44 x 10^^ free protons for 20 kt. The differential flux of 
reactor neutrino d(f){E)/dE can be computed as. 


d(t){E) 

dE 




(B.2) 


where Pth is the thermal power of the reactor, {E) ~ 210 MeV is the average energy 
released by per fission computed by taking into account the ratios of the fuel compositions 
of the reactor (see below). 


We can replace, in a good approximation, the reactor complex consisting of 6 and 4 
reactors, respectively, at Yangjiang and Taishan sites by a single reactor with the thermal 
power of 35.8 GW placed at the baseline L = 52.5 km from the JUNO detector. We also 
include the contributions from the far reactor complexes at Daya Bay (with the baseline of 
215 km) and Huizhou (with the baseline of 265 km) sites, which contribute, respectively, 
about 3% and 2% in terms of the total number of events. 


For the reactor spectra S{E) we use the convenient analytic expressions found in [38] 
with the typical fuel compositions of the reactors, ^^^U: ^^®Pu: ^^®U: ^^^Pu = 0.59: 0.28: 
0.07: 0.06, obtained by taking time average of Fig. 21 of [48]. We note that S{E) is nothing 
but the number of neutrinos being emitted per fission per energy (MeV). 

Furthermore, we also include, for completeness, geoneutrinos coming from the decays of 
U and Th inside the earth in the same way as done in [31], despite that it is not important for 
our main purpose. Assuming the input (true) geoneutrinos fluxes of (/>(U) = 4.0 x 10® cm“^ 
s“^ and 0(Th) = 3.7 x 10® cm“^ s“^, the expected U and Th geoneutrinos induced events 
are, respectively, ~ 1.9 x 10^ and ~ 5.4 x 10^ for 5 years of exposure at JUNO detector. 
For simplicity, we consider only the averaged standard oscillation for geoneutrinos, and 
ignore the decay effect for them. Since the presence of geoneutrinos has minor impact on 
the standard oscillation study at JUNO, their decay effect is even minor and we believe 
that this is a fairly good approximation. In fact we have verified that the presence of 
geoneutrinos has virtually no impact on the sensitivity to the decay effect. In the fit, 
we let the fluxes of geoneutrinos vary freely subject to the pull terms in Eq. (4.5) with 

(Tgy — — 20%. 

R{Ee, Evis) is the function which takes into account the finite energy resolution of the 
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detector and is given by 


R{Ee, E^is) 


1 

,_ _exp 

V^a{Ee) 


1 / -j- TTTg E^i] 

2 V cr{Ee) 


2 


(B.3) 


where the energy resolution is assumed to be [49], 

(^{Ee) ^ 3% (1 + rj) 

{Ee + nie) y/{Ee + me)/MeV' ^ 

where rj is introduced to take into account the uncertainty in the energy resolution. The 
energy resolution in (B.4) amount to consider only the stochastic term. 

The expected total number of events at JUNO for the 5 years of exposure with 100% de¬ 
tection efficiency is 1.41 x 10^. The individual contributions from medium-baseline reactor 
sites (at Yangjiang and Taishan), far reactor sites (at Daya Bay and Huizhou), geoneutrinos 
are, respectively, 1.39x10^, 6.70x10^ and 2.40x10^. 


C Correlations among the oscillation, systematic and decay parameters 

Here, we discuss the correlations among the mass and the mixing parameters, system¬ 
atic uncertainty parameters, and the decay parameter Tsjm^. For this purpose we have 
examined all possible correlations among the mixing parameters and the uncertainty pa¬ 
rameters.^^ The calculation is done under the same conditions (running time etc.) as 
assumed for calculating the lifetime bound given in Fig. 2. 

The general features of the correlations among the oscillation parameters, systematic 
uncertainty parameters and the decay parameter (ms/rs) can be summarized as follows: 

• The effect of decay is visible only in the contours which involve the energy resolu¬ 
tion (T£;(l -|- ry), 013 and Am|i, as implied by the results shown in Table 1.^^ 

• The correlations between r^jms and aEi^+r]), 0i3 exist, but in easily understandable 
way. 

C.l Decay lifetime vs. 0i3 and the energy resolution 

In this subsection, we highlight only the correlations between the 1/3 lifetime and (a) sin^ 0i3 
and (b) the uncertainty in energy resolution. All the other combinations not shown here 
and in section 6.2 will be shown in the next subsection, see Fig. 6. 

They are the following 32 combinations (to be shown in Fig. 6 ) in addition to those already shown in 
Fig. 4 and the one to be shown in Fig. 5: 

(a) sin^ 6 ' 12 -sin^ @ 13 , (b) sin^ @i 2 -Am 3 i, (c) sin^ @i 3 -Am|i, (d) Amli-Amli, 

(e) sin^ @12-(1 -I- Creac), (f) Sin^ @ 13 -(1 -f ^reac ) , (g) Amii-(1 -f ^reac), (h) Amli-( 1 + ^reac), 

(i) sin^ @ 12-(1 -I- Cu), (j) sin^ @i 3 -(l + ?u), (k) Amii-(1 -|- ^u), (1) Am 3 i -(1 -|- ^u), 

(m) sin^ @i 2 -(l + ^Th), (n) sin^ @i3-(l + Cxh), (o) Amii-(1 + ^xh), (p) Amii-(1 + ^xh), 

(q) sin^@i 2 -crB(l+ ? 7 ), (r) sin^ @i 3 -ctb( 1 -|- 77 ), (s) Amii-(TB( 1 -|- 77 ), (t) Am 3 i-cr£;(l-|- 77 ), 

(u) (1 -t- ^u)“(l “t“ ^reac), (v) (1 + ^Xh )-(1 “t“ ^reac) (w) (1 ^u)“(l “t ^Th), (x) 01^(1 -f 77)-(l -f ^reac) 

(y) (1 + ^u)- 0 £;(l-b r;), (z) (1-b CTh)-o-£;(l-b 77 ) (A) sin^ @i 2 -r 3 /m 3 , (B) Am^i-Tilms, 

(C) Amji-Ta/ms, (D) (l-b^ reac )-T3/m3, (E) (1-b Cu)-'r 3 /m 3 and (F) (1-b CTh)- 7 - 3 /m 3 . 

If we remove all the priors for the mixing parameters (i.e., pull terms) the feature changes dramatically. 
There arises an extended allowed region up to sin^ @13 ~ 0.044 at 2a C.L. with a larger flux normalization 
by about 5%. 
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Figure 5. The regions allowed at 1, 2, and So CL are drawn, respectively, by the filled colors 
of blue, yellow and green in the parameter space spanned by (a) sin^ 613 - T 3 / m 3 (left panel) and 
(b) aEi^ + rj) ~ (right panel) for the case where no decay is assumed in input. The case 

where Tslrris = s/eV is assumed in input is also shown by the black dotted curves. 5 years 

of exposure is assumed. 


In Fig. 5, we show the correlations in space spanned by (a) sin^ ^13 - T^jm^ (left panel) 
and (b) (Te(1 + ??) - Tsfms (right panel). We notice that there are correlations, i.e., the 
effect of decay, both in the left and right panels, but only in near the lower end of the 
allowed region of the lifetime. From the left panel of Fig. 5 we learn that the decay effect 
(diminishing the amplitude of the wiggles) can be compensated to some extent by making 
013 larger, which produces a new allowed region toward the large 0i3 direction. 

Corresponding to the newly emerged region, there also arises a new allowed region in 
the iT£;(1 + rj) - Tzjraj, space, as seen in the right panel of Fig. 5. It is located at near the 
top of the peninsula, a slightly distorted region toward small (Te (1 + r\) direction in the 
right panel of Fig. 5. The shape reflect a better energy resolution in the newly allowed 
region. We note that due to the above sin^ 0i3 - crE(l + ?/) correlation, there always exist a 
slightly extended allowed region in any one of the correlation plots which involve sin^ 613 
or (7^(1 + T]) toward the direction of larger 0i3 and smaller (T£;(1 + r/). 

C.2 Correlations of parameters: miscellaneous cases 

We now show, for completeness, the allowed regions spanned by all the remaining 32 
combinations of parameters considered in this work (which are described explicitly in the 
footnote 12), except for those already shown in Figs. 4 and 5. The meanings of the filled 
colors and line types are the same as in Figs. 4 and 5. 

As we can see from these plots that, in general, there are no significant differences 
among the allowed regions shown by the filled colors (no decay in input but allowed in 
the fit), solid black curve (standard oscillation fit without decay) and black dashed curves 
(decay effect both in the input and in the fit). Therefore, as a whole, the impact of the 
decay is rather small. Furthermore, we note that there is no significant newly induced 
correlations due to the decay effect. 
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Figure 6. Allowed regions in the space spanned by all the combinations of parameters not shown 
in Sec. 6 and in the previous subsection. The filled color regions correspond to the case where no 
decay is considered for input but allowed in the fit. The black solid curves correspond to the case 
of the standard oscillation fit without decay whereas the black dotted curves correspond to the case 
where the true value of Ta/ms is assumed to be 10“^° s/eV. 
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